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An improved one-dimensional mathematical model based on Pulsed Flow Equations 
(PFE) is derived by integrating the axial component of the momentum equation over the 
transient Womersley velocity profile, providing a dynamic momentum equation whose co- 
efficients are smoothly varying functions of the spatial variable. The resulting momentum 
equation along with the continuity equation and pressure-area relation form our reduced- 
order model for physiological fluid flows in one dimension, and are aimed at providing 
accurate and fast-to-compute global models for physiological systems represented as net- 
works of quasi one-dimensional fluid flows. The consequent nonlinear coupled system of 
equations is solved by the Lax-Wendroff scheme and is then applied to an open model 
arterial network of the human vascular system containing the largest fifty-five arter- 
ies. The proposed model with functional coefficients is compared with current classical 
one-dimensional theories which assume steady state Hagcn-Poiscuille velocity profiles, 
either parabolic or plug-like, throughout the whole arterial tree. The effects of the non- 
linear term in the momentum equation and different strategies for bifurcation points 
in the network, as well as the various lumped parameter outflow boundary conditions 
for distal terminal points are also analyzed. The results show that the proposed model 
can be used as an efficient tool for investigating the dynamics of reduced-order models 
of flows in physiological systems and would, in particular, be a good candidate for the 
one-dimensional, system-level component of geometric multiscale models of physiological 
systems. 
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1. Introduction 

Extensive research has been done on modeling human physiology. Most of the work 
has been aimed at developing detailed models of specific components of physiological 
systems such as a cell, vein, molecule, or heart valve. One example of such an effort 
is the three-dimensional (3D) computational model of the heart developed at NYU 
by Peskin and McQueen^. The challenges of the interaction of mechanics with 
biology and medicine have recently been outlined by Bugliarello- 2 . In particular, 
modeling of internal flows through elastic vessels has also been studied intensively 
over yearsP"^, and a brief history of the arterial fluid mechanics can be found in 
a recent paper by Parker -. Wide range of the proposed models may be classified 
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by their dimensionality. As far as those models are concerned, depending on the 
phenomena to be studied, numerous degrees of simplifications can be considered 
going from fully 3D fluid-structure interaction models'sH^ to the simplified reduced- 
order models 12 " 14 !, and from complex non-Newtonian fluid modelg 15 " 20 ! to idealized 
Newtonian fluid models^^l. 

There has been recent interest in large-scale fully 3D computations of patient- 
specific human arterial tree because of their importance in the emerging field of 
personalized medicine. This 3D numerical analysis, which is generally based on 
MM imaging of the patient's vasculature, provides a very detailed view of local 
flow features such as recirculation and shear stress in a patient's circulatory system, 
but is enormously computationally intensive. Grinberg et al. report in^that simu- 
lating one cardiac cycle in an arterial tree that included the largest arteries required 
27.7 hours of computational time using 40,000 processors. While fully 3D simula- 
tions such as these are invaluable to our understanding of human physiology, and 
to treating common circulatory system disorders such as arterial stenosis and heart 
valve malfunction, most hospitals and clinics do not have ready access to supercom- 
puting clusters, and an approach that is still very accurate, but less computationally 
intensive is desired in the clinical or hospital setting. Global computational models 
of patients' entire cardiac functioning that can give useful results in a reasonable 
amount of time on a desktop computer are needed. Some of geometric multiscale 
models of the cardiovascular system have been recently developed that join fully 
detailed 3D simulations of a region of interest in a patient, for example, the carotid 
artery for a carotid artery stenosis, with global, macroscopic level reduced-order 
models for the rest of the patient's cardiovascular system, thereby reducing the 
spatial and computational complexity of the modeP^IHSZJ 

The simplest reduced-order models for patient vasculature are spatially homo- 
geneous, zeroth-order (0D) lumped parameter models based on an analogy between 
the cardiovascular system and an electrical circuit. In this approach, each segment 
of the system is modeled by basic circuit elements yielding a set of ordinary dif- 
ferential equations in the time domain. The complex network that comprises the 
human circulatory system can be simulated efficiently in this The main 

disadvantage of these 0D models is that they are unable to take into account some 
important features of cardiovascular function, such as wave propagation along ves- 
sels, which is patient-specific and critical to the accuracy of the detailed, microscopic 
level models. 

Another approach to reduced-order modeling of the cardiovascular system is 
to develop spatially one-dimensional (ID) linear or nonlinear models consisting 
of coupled partial differential equations that describe the behavior of the venous 
cross sectional area, the averaged blood velocity, and the averaged blood pressure 
throughout the system^ 6 * 12 ! 31 * ^. These reduced-order ID modeling efforts offers an 
excellent compromise with anatomic accuracy and has a great potential in clinical 
and research applications?, as well as in the other mechanical systems such as 
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valveless impedance pumping for open or closed configurations^. They are based on 
the assumption that the flow velocity along the elastic vessels is much greater than 
the flow velocity perpendicular to the longitudinal axis, resulting a nonlinear coupled 
system of averaged continuity and momentum equations. An additional equation is 
used to close the system by providing the constitutive relationship between the 
cross-sectional area at a point along the system, and the pressure at that point. In 
the geometric multiscale approach to modeling the cardiovascular system, these ID 
models of a patient's vasculature are often truncated after a number of arterial tree 
generations, and supplemented with appropriate boundary conditions provided by 
0D lumped parameter models that model the effects of the rest (the distal portion) 
of the vasculatures™™ 

As the ID modeling and computations have become more or less routine, atten- 
tion has turned to improvement the basic theory, especially in the context of bifurca- 
tions^!, elastic wall modeling * 47 * 48 !, and permanent constrictions (i.e, stenoses^) or 
forced constrictions (i.e., valveless pumping^). In the ID theory, momentum equa- 
tion is derived by integrating an assumed velocity profile over the cross-sectional 
area of the tube. In all ID modeling studies cited above, the assumed velocity pro- 
file has either plug-like form (corresponding to large Womersley numbers for larger 
vessels) or parabolic form (corresponding to small Womersley numbers for smaller 
vessels). In this classical formulation, the effect of changing shape/frequency param- 
eter (Womersley numbers) on pulsatile flow is neglected due to the approximated 
steady state Hagen-Poiseuille velocity profile giving a momentum equation with 
constant coefficients. In this study, a new one-dimensional model is derived by inte- 
grating the momentum equation over the transient Womersley flow velocity profile. 
Therefore, the coefficients in the averaged momentum equation changes continu- 
ously via Womersley numbers allowing us to simulate global physiological systems 
more appropriately. 

Many numerical methods have been proposed to solve hyperbolic system of 
equations can be used for ID modeling, for example, method of characteristics^!, 
finite differenc e ^ 34 * 35 * and finite/spectral element method^ 33 * 42 !. The two-step Lax- 
Wendroff scheme in the finite difference methods and its finite element equivalent 
form, Taylor-Galerkin scheme are common examples of second-order accurate stan- 
dard algorithms for computing these flows. The finite-difference discretization with 
Lax-Wendroff scheme is used in this study to solve the improved reduced-order 
pulsed flow equations. It is second-order accurate in time and space. The PFE can 
accommodate multiple branching, nonlinear elastic wall effects, tapering vessels, 
some two-dimensional viscous and inertial effects, and different boundary condi- 
tions, as well as the effects of external forces such as gravity and outperforms the 
classical one-dimensional theory for frequecy/shape parameter changing internal 
flows such as flows in arterial tree. 

The outline of this paper is as follows. In section 2, we first introduce the Pulsed 
Flow Equations (PFE) by extending the classical ID theory with assumed Worn- 
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ersley velocity profiles changing continuously via a shape parameter. The different 
strategies for boundary conditions at terminal and bifurcation points are also pre- 
sented in this section. The Lax-Wendroff scheme is introduced in section 3 in order 
to solve PFE efficiently. In section 4, the model and algorithm is then applied to 
an open model of a typical arterial network containing the largest fifty-five arteries 
with proximal boundary conditions at the Ascending Aorta, coupled with a lumped 
parameter model at the distal points of the vasculature to provide outflow boundary 
conditions. The relative importance of modeling strategies such as the influence of 
the assumed velocity profiles in ID modeling, various models for bifurcation and 
terminal point boundary conditions is also analyzed and compared each others. 
Finally, in section 5, we present some conclusions and perspectives related to this 
work. 



2. Mathematical Modeling 

2.1. Derivation of the Pulsed Flow Equations 

The mathematical model is based on the mass and momentum conservation equa- 
tions of incompressible fluids streaming through the tubes, assuming axisymmetric 
flow of constant fluid density, temperature, and viscosity. The continuity and axial 
momentum equations of incompressible Newtonian fluid (a valid assumptions for 
blood flows through the large and medium vessels'^) are given in axisymmetric 
cylindrical coordinates as: 

d(rv) d(ru) 



dr dx 



(1) 



du du du I dp . 1 d , du. d 2 u. 

— +u— +v— = — jf +v (- — r— + -j—) 2) 

at dx or p ox r Or Or dx z 

where r and x are the coordinates in radial and axial directions. Here, p is the 
density, v is the kinematic viscosity, p is the pressure, u and v are the velocity 
components in axial and radial directions. Since the blood density is constant and 
the flow in confined in the artery, the gravity effect can be included in the pressure 
term. We seek to derive a reduced-order system of equations that govern the aver- 
aged axial velocity, u, the pressure, p, and the cross-sectional area, A, in internal 
biological fluid flows, while capturing some essential effects occurring in biologi- 
cal fluid flows that are not captured by the standard equations for incompressible 
fluid flows in one-dimension. This can be accomplished by integrating the equations 
over the cross-sectional area. Integrating continuity equation over the cross-section 
results in: 

OA d(uA) 



dt dx 



(3) 



where we used v(R) = Here the cross-sectional area, A = ttR 2 , and the averaged 
velocity, u = J u(r, x,t)rdr, are functions of one-dimensional space variable x, 
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and time variable t. In order to integrate the axial momentum equation over the 
cross-section, an assumed velocity profile need to be taken. It is assumed that the 
velocity profile has the form: 

u(r,x,t)=f(r)u(x,t) (4) 

In the classical one-dimensional models, f(r) is taken as steady state Hagen- 
Poiscuillc (HP) profile: 

/HP (r) = 7 + 2 r 
7 R 

where the coefficient 7 = 2 corresponds to fully developed, steady, incompressible, 
Newtonian flow in a cylindrical tube, which is also called Poiseuille flow, and is used 
in many studiest^MEI!]. Measurements show that the parabolic velocity profile is 
not quite accurate for the larger arteries. Alternatively, the plug-like velocity profile 
with 7 = 9 is extensively applied for the arterial tree networlp2Hl!l. Here, instead 
of these approximations, we use the following profile corresponding to pulsatile 
Womersley flow^Sl 

where 3? represents the real part of the complex variable and Jo, and J\ are the 
zeroth and first order Bessel functions of first kind^. Here A is defined as: 

A=(^)Wo, Wo=^Ji?(*,t) (7) 

where i 2 = — 1, uj is the angular frequency of transient flow and Wo is Womersley 
number which can be thought as frequency or shape parameter (representing the 
various size of tube radius). Integrating Eq. [2] over the cross-section by using Eq. [4] 
with Eq. [6j the averaged momentum equation becomes: 

du ^ _du 1 dp ^ d 2 u u 
dt dx p dx dx 2 A 
where the coefficients becomes: 



where the asterisk represents the complex conjugate. In above formulation, we have 
used the identities: 

dJ (r) , - drJi(r) 

-fir- = -Mr), ~^-=rMr) (11) 
J (^A)rdr = ^ J X (A) (12) 
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J*j2(L A)rdr= ^ { j2 {A) + j2 {A)] (13) 

J*J l (LA)dr=*[l- MA)] (14) 

On the other hand, by assuming Hagen-Poiseuille velocity profiles the the coefficients 
of momentum equation, Eq. [8j are obtained as: 

a= I ±|, /3 = 2( 7 + 2) (15) 
7 + 1 

which are constants either 7 = 2 for parabolic velocity profiles, or 7 = 9 plug-like 
velocity profiles. However, in our formulation given in Eqs. [9 |l0[ the coefficients of 
momentum equation become functions of Womersley number providing a continuous 
dynamical partial differential equation for flows through tubes with various sizes in 
diameter as occurred in the physiological systems. For example, corresponding Wo 
numbers for some arteries from largest to medium sizes in radius is shown in Table 
1 where the blood is considered to have a constant dynamic viscosity of 0.0045 
Ns/m 2 , and a constant density of 1050 kg/m 3 with having 1 s oscillatory cycle 
period corresponding to 60 heart beats per minute. 



Table 1. Womersley numbers for some vessels in the arterial tree network. 



Arteries 


L 


Rtop 


Rbot 


CO 


Wo 




(cm) 


(cm) 


(cm) 


(m/s) 


IT = Is) 


Ascending Aorta (1) 


4.0 


1.470 


1.440 


4.30 


17.8 


Thoracic Aorta I (18) 


5.2 


0.999 


0.999 


4.30 


12.1 


Abdominal Aorta I (28) 


5.3 


0.610 


0.610 


5.00 


7.38 


Superior Mesenteric (34) 


5.9 


0.435 


0.435 


5.10 


5.27 


L. Carotid Artery (15) 


20.8 


0.370 


0.370 


5.30 


4.48 


R. Femoral Artery (52) 


44.3 


0.259 


0.190 


8.60 


3.14 


Celiac B (30) 


1.0 


0.200 


0.200 


7.30 


2.42 


L. Interosseous (24) 


7.9 


0.091 


0.091 


14.3 


1.10 



The normalized Womersley velocity profiles u(r,x,t) /u(x,i) = f(r) are plot- 
ted in Fig. [T] The Hagen-Poiseuille velocity profiles are also shown in the figure. 
For smaller Wo number, (Wo < 1), there is no significant difference between the 
Womersley and Hagen-Poiseuille parabolic profile. The Eq. [9] and Eq. [10] for coeffi- 
cients of nonlinear and viscous drag terms in the momentum equation with varying 
Womersley number can be also reduced to the coefficients for classical ID theory 
(7 = 2) for small Womersley numbers. For convenience if we define, e = Wo/A 
then A = e(i— 1). By using the asymptotic formulas for Bessel functions with small 
arguments^, we have 

Jo(A)~l + *^ (16) 
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Fig. 1. Steady Hagen-Poiseuillc parabolic (7 = 2) and plug-like(7 = 9) velocity profiles and five 
representative snapshots of the peak velocity profiles of the transient Womersley flow. 




Fig. 2. Varying coefficients of (a) nonlinear term and (b) viscous drag term in the momentum 
equation. 



■MA)~f[-(l+^) + i(l-^)] 



(17) 
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Substituting Eqs. l"6ff"7 into the Eqs.[9p0l we get the coefficients of classical theory: 



/3~8 (18) 

Since in the arterial tree Womersley number is bigger than that, therefore the classi- 
cal theory has been extended for the plug-like velocity profile (7 = 9) as mentioned 
before. This assumption approximates the velocity profile better for larger arteries, 
but takes one constant value all larger vessels. We improved one-dimensional theory 
by using the Womersley velocity profiles with functional coefficients, a and j3 given 
in Eqs. [9] |10[ of momentum equation which are shown in Fig. [2j Since Womersley 
parameter is changing through the tube segments in the network, the proposed 
functional form of a and j3 is changing according to the Eqs. [9]|T0| through the 
network. This improvement in the one-dimensional theory allows us more accurate 
simulation of network type flows as occurred in human arterial tree for varying ves- 
sel radius. As mentioned earlier, previous studies for flows in human arterial tree 
have been performed using constant values for a and /3 through the whole network. 
The influence of velocity profiles in ID modeling is analyzed in section 4. 



2.2. Constitutive relationship 

In order to close the system of equations, the constitutive relation between the 
pressure and the cross-sectional area must be specified at a point in the system. 
Since the viscoelastic effects are small within the physiological range of pressure, 
the linear theory of elasticity can be used* 6 * 12 * 29 * 47 !. According to Hooke's Law for 
thin walled tubes, the tangential and axial stresses and strains are related through 
the following formulas^! 

s t = % - ^ (19) 
EE K ' 

&s vat 

£ °=E-Jf ™ 

where E is the Young's modulus of tube material and v is the Poisson ratio, which 
is equal to 0.5 for incompressible materials. Considering that the elongation of the 
tube is zero (e s = 0), then the tangential strain becomes: 

For a thin walled tube, the tangential tensile stress acting on it is 

vt = j(p-Po) (22) 

where po is the system's reference equilibrium pressure {p is the fluid pressure at 
inner surface, and po is the pressure at outer surface), and h is the tube wall thick- 
ness. By using the definition of tangential strain, e t = (2ttR — 2n R ) / '(2n R ) , and 
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comparing it with Eq. (21), an algebraic expression of the pressure versus tube 




cross-sectional area is obtained: 

p = Po + ^n ^(W^) ( 23 ) 

R [l-v 2 ) 

where Aq is the reference equilibrium area of the tube (the area when the p = po)- 
The Young's modulus of the thin-walled tube is related the wave propagation speed 
by the Moens-Korteweg equation^ 



/ Eh 

00 = V w^j (24) 

where cq is the wave propagation speed which depends on the properties of the 
elastic tube. Then, the constitutive relationship between pressure and area can be 
rewritten in the form of: 




p=p + 2p{c Y(l-\ /-H (25) 



As shown in Table 1, the corresponding wave propagation speed is in the range 
from 4 to 15 m/s and increases by decreasing vessel radius due to higher Young's 
modulus of smaller vessels. To account tapering effect for some of the vessels, the re- 
lationship for the radius along an individual artery is often specified in the following 
exponential relationship! 13 * 29 ! : 

_/„ I R bottom \ X 

Ro(x)=R top e tn( ^ )l (26) 

where Rtop and Rbottom are the reference proximal and distal radii and L is the 
length of the individual segment. The proximal (top) and distal (bottom) radius 
values under some reference condition for each particular segments for the human 
arterial system are well-documented in the literatur e ! 13 * 34 * 37 * 4 1 ^ In this study, the 
reference equilibrium radius values along the tubes through the network are com- 
puted by using Eq. |26j Then the local Womersley number and the corresponding a 
and /3 coefficients through the whole network can be pre-computed. 



2.3. Boundary conditions 

In order to finalize mathematical model, an appropriate boundary and initial con- 
ditions should be specified to the problem. For pulsatile flow problems, initial con- 
ditions does not have effect on the final solution. The results should convergence a 
periodic state after a few dummy simulation of oscillation cycles. Usually the hy- 
perbolic system of equations has a positive wave propagation velocity much greater 
than the actual velocity of flow (cq ^ u). This is also true for blood flow compu- 
tations in the arterial network. Consequently, we only need to specify either the 
area (pressure) or velocity values at the boundaries. Since the area and pressure are 
related each other by an algebraic constitutive relationship, Eq. [25j it is enough to 



December 4, 2012 2:37 WSPC/WS-JMMB ws-jmmb 



10 Omer San and Anne E. Staples 

specify pressure boundary conditions to the problem. For the open system of net- 
work flows, as studied in many arterial tree computations, the following boundary 
conditions must be established: (i) An equation at the inlet to the network; (ii) an 
equation at the outflow from each of terminal points of the network; (iii) a model 
at the each of the bifurcations. 



2.3.1. Inflow boundary condition 

As an inflow boundary condition, we enforce either a time dependent area function, 
Ai n (t), or a time dependent velocity, Ui n (t). In most blood flow simulations the 
flow rate, Qi n (t) = uA, or pressure, p% n (t) waveform can be specified as an inflow 
boundary condition at the proximal end of the network (i.e., the point where the 
arterial system is closest to the heart). In this work, we consider that the flow rate, 
Qi n (t) , is given as inflow boundary condition. The other variables can be computed 
with the help of discrete version continuity equation: 

= A » ~ Ix (<+1j4 " +1 ~ (2?) 
where subscripts refer to spatial index, and superscripts denotes the discrete time 
level. The subscript locates the inlet point and the subscript 1 locates the nearest 
discrete point to the inlet. Then, the corresponding velocity condition becomes 
trivial (v,i n = Qi n /^o + ) an< ^ the i me t pressure is determined from Eq. 
used this methodology for all the cases presented here. 



25 



We have 



2.3.2. Outflow boundary condition 

For an open system arterial flow, the model is usually terminated using a lumped 
parameter impedance model for practical reasons. In the geometric multiscale ap- 
proach to modeling the cardiovascular system, the ID models of a patient's vascu- 
lature are often truncated after a number of arterial tree generations, and supple- 
mented with appropriate boundary conditions provided by 0D lumped parameter 
models that model the effects of the rest (the distal portion) of the vasculature. The 
well accepted three-element windkessel model shown in Fig. [3|a) with two resistors 
and a capacitor is used at the terminal points^ 30 * 34 l In the windkessel approach, 
the resistor elements represent the hydraulic resistance of the tube and the conduc- 
tor element represents the elastic capacity of the tube. Using the circuit law, the 
corresponding relationship between the pressure and flow is given by: 

d(uA) 1 dp | p B.! uA 

dt R\ dt R\R.2Ct R2 R\Ct 

where R\ + i?2 = Rt is the total resistance and Ct is the total compliance of the 
terminal branches. The major advantage of this model is that it accounts for both 
resistant and compliance effects of the distal vessels beyond the point of termina- 
tion. Alternatively, a simpler pure resistance model shown in Fig. [3jb) that does 
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Fig. 3. Terminal (outflow) boundary conditions: (a) the three-element windkessel model with two 
resistor and one conductor, (b) pure resistor model. 



not consider compliance effects beyond the terminal points can also be used as an 
outflow boundary condition with the corresponding relationship given as: 

a = 56 < 29 > 

These two lumped-parameter models are applied as outflow boundary conditions 
for terminal points and compared with each other in this study. 

2.3.3. Bifurcation point boundary conditions 

The one-dimensional theory derived for a flexible tube can be extended to handle 
the network flows by imposing a suitable interface condition at the branching points 
of the network. Assuming that all bifurcations occur at a point with one mother 
segment and two daughter segments as shown in Fig. |4j the following nine quan- 
tities, A (1 \ p (1 \ u (2) , A (2) , p (2) , m (3) , A( 3 \ p (3) , need to be specified at each 
bifurcation points. Since the pressure values are algebraically related to the area 
values via Eq. |25[ there are only six independent unknown quantities for each bifur- 
cation point in the network. The superscripts here indicate the segment numbers: 
(1) represents the mother tube, and (2) and (3) denotes the daughter tubes. Flow 
conditions at the bifurcation points can be represented by conservation of mass and 
pressure in the following formal 

U a) A (l) =^(2) j4 (2) +iZ (3) A (3) (30) 

pW=pW (31) 

pW=p& (32) 

Alternatively, the Eq. [31] and Eq. [32] for pressure equalities can be improved by 
using the total pressure equalities as 3 ^* 

+ \p{u^f =p< 2 > + \p(u^r (33) 
P^ + \p(u^f=p^ + \p(u^f (34) 
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Fig. 4. The schematics for the bifurcation boundary conditions. Superscript 1 denotes the mother 
segment and superscripts 2 and 3 denote the two daughter segments. The subscript * represents 
the values for the nearest solution points to the boundary points. 



The equations above provide three equations for six unknowns for each bifurcating 
point. Three additional equations come from either Riemann invariant d 33 * 35 ! or dis- 
crete continuity equations^. Looking at the problem from a characteristics point of 
view (using Riemann invariants), information can only reach the bifurcation point 
from within mother tube by a forward-traveling wave. Similarly, within the daughter 
segments, information can only reach the bifurcation point by a backward-traveling 
wave. For PFE equations, Riemann invariants are approximated as: 

q ± =u±±c [l-(A /A) 1 / 4 ] (35) 





+ 4c«[l- 


(4 X) M (1) ) 1/4 ] 


= qt 




-44 2 >[1- 


(4 2) M (2) ) 1/4 ] 


= q* 




-44 3) [i- 


(A^ /A {3) ) 1/4 } 


= q* 



The derivation of Eq. 35 is shown in Appendix. Using the positive Riemann invariant 
for the mother tube and the negative invariants for the daughter tubes, the three 
additional equations are: 

(36) 
(37) 
(38) 

where are values that can be easily computed at the nearest solution points in 
the discretized system of equations. 

Alternatively, the discrete version of the continuity equations can be used as 
closure equations. For the mother segment, it is given as: 

— + A * * = (39) 

At Ax y ' 

Similarly for the daughters: 

4(2) _ A (2)n u i2) A {2) - u^AW 

+ u * A * u = o (40) 

At Ax 
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(41) 



At Ax 
where the subscript * represents the nearest discrete point to the bifurcation point 
in the vessel. The superscript n represents the previous time-step values. With com- 
bination of the strategies described above we can construct four different treatments 
for the bifurcation points. The summary is given in Table 2, and the effect of these 
procedures to the global model will be discussed in Section 4. The six unknowns 
are then solved with the Newton-Raphson algorithm^ using a nonlinear set of six 
equations. We use the values from the previous time step as an initial guess in the 
Newton-Raphson algorithm. The Newton-Raphson algorithm converges quadrati- 
cally in a few iterations, usually one or two iterations per time step. 



Table 2. Strategics for bifurcation boundary points for computing six unknowns. 



Procedure set 


Set of equations 
















Characteristics & total pressure 
Discrcatc continuity &; total pressure 
Characteristics & pressure 
Discreate continuity & pressure 


Eq. 
Eq. 
Eq. 
Eq. 


3(5 
3!) 
36 
39 


Eq. 
Eq. 
Eq. 
Eq. 


37 
TO 
37 
TO 


Eq. 
Eq. 
Eq. 
Eq. 


38 
11 
38 
11 


Eq. 
Eq. 
Eq. 
Eq. 


30 
30 
30 
30 


Eq. 
Eq. 
Eq. 
Eq. 


31 
31 
33 
33 


Eq. 
Eq. 
Eq. 
Eq. 


32 
32 
31 
31 



3. Lax-Wendroff Scheme 

As stated earlier, the PFE can be solved by classical numerical methods for hyper- 
bolic partial differential equations such as the MacCormack or Lax-Wendroff finite 
difference schemes or Taylor-Galerkin formulation in the finite element framework. 
Governing PFE that describe the temporal and spatial evolution of pressure, ve- 
locity and cross-sectional area along the elastic tubes and can be rewritten in the 
following vector form: 

dQ dF_ 

dx 



where 



'A' 






,F = 


u 



dt 
Au 

fw 2 + 2( Co ) 2 (l- 



= S 



(42) 



,S = 



1 dpo 
p dx 



dx 2 



The problem given in Eq. [42] can be solved by Lax-Wendroff scheme using the 
following two-step procedure^ 



Q 



i+i 
n+1 



Qi = Q\ 



(43) 
(44) 
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where superscripts and subscripts show the time index and spatial index, respec- 
tively. Inside each tube, the n+1 time-step values are obtained from the n time-step 
values using two subsequent steps. The inflow/outflow boundary conditions and the 
internal boundary conditions at the bifurcation points are provided in the previ- 
ous section. This scheme is second-order accurate in space and time. The stability 
criterion for the Lax-Wendroff scheme is 



where cq is the speed of wave propagation defined earlier by the Moens-Korteweg 
relationship. 

4. Results for arterial network 

In this section a representative model for the human cardiovascular system con- 
taining the largest fifty-five arteries is simulated. The schematic description of this 
network flow is illustrated in Fig. [5j The physiological data including length, proxi- 
mal and distal radius of every fifty-five segments and total resistance and compliance 
coefficients of terminal points and wave propagation speed (including the elastic pa- 
rameters of the arteries) were prescribed based on the data reported j n IM3zEy] The 
blood is considered to have a constant dynamic viscosity of 0.0045 Ns/m 2 , and a 
constant density of 1050 kg/m 3 . The equilibrium pressure is chosen as Po = 98 
mmHg, which is assumed as the reference mean pressure corresponding to the equi- 
librium reference diameters of the vessels^. The flow waveform is specified at the 
proximal end of the network (at the beginning of the Ascending Aorta) by using 
the first ten Fourier harmonics^ with the period of one second, which is plotted 
in Fig. [6] corresponding to 60 heart beats per minute giving approximately 5.18 
mL/min cardiac output to the arterial tree. 

Although this is a time dependent problem, exact initial conditions can not be 
imposed at all discrete points except at the proximal and distal boundary points 
explained in Section 2. Fortunately, the viscous effects damped out the discrepancies 
due to incorrect initial pressure and flow distribution, and in case of periodic flow 
the solution converges within two or three cycles as shown in the Fig. [7] The flow 
waveform is periodic (specified Fourier harmonics), and the pressure waveform at 
the inlet is converging after the third cycle. This convergence is true for all the 
segments in the network. Therefore we will present our converged results in rest of 
the paper. 

4.1. Grid dependance study 

We first investigate the effect of the the spatial resolution within arterial tree net- 
work with a small time step chosen as At = 2.5 x 10~ 5 s. Comparisons of different 
resolutions on pressure and flow waveforms are shown in Fig. [8] and Fig. [9] for the 
distal point of the Abdominal Aorta I (segment # 28), which is typical of a large 



CFL = 



At(u + c ) 
Ax 



< 1 



(45) 
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Fig. 5. The arterial tree model including the largest fifty-five arteries. 




Fig. 6. The periodic inflow for the proximal end of the Ascending Aorta (segment # 1). 

artery, and for the distal point of the R. Femoral (segment # 52), which is repre- 
sentative of a typical smaller sized artery. In figure labels, N represents the total 
number of grid points in the network. We can conclude that the Ax = 2 mm along 
the arterial tree is enough to capture all the waveform correctly. This was also 
outlined inP^l The computed waveforms are similar to those computed in^, and 
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Fig. 7. Convergence of pressure and flow waveforms at the inlet. The thin blue line shows the 
specified periodic inflow boundary condition, and the thick red line shows the pressure waveform. 



exhibit the general characteristics of the systemic circulations observed in experi- 
mentally^. The shapes of the pressure and flow waveforms were generally in good 
agrement with experimental measurement and the effects of the other mathematical 
modeling issues will be analyzed systemically in the following subsections. 
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Fig. 8. Effects of varying grid resolutions on the pressure waveform. 
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Fig. 9. Effects of varying grid resolutions on the flow waveform. 



4.2. Effects of the nonlinear term in the momentum equation 

In order to quantify the effect of nonlinear term in the momentum equation, we 
perform two numerical runs such that one with a given in Eq. [9]and one with a = 0. 
Linear model assumes a — 0, and for nonlinear model, the a is changing through the 
arteries in the network and it is determined from Eq. [9j For both model, the other 
conditions and parameters are the same, i.e., the viscous drag term coefficient j3 is 
computed from Eq. [10] Computed pressure and flow waveforms are given in Fig. [10] 
and Fig. [TTJ respectively. The model with nonlinear term estimates a little bit bigger 
pressure amplitude through the arterial tree. Since the wave propagation speed is 
much bigger than the actual velocity this small difference is expected. It can also 
be seen that there is very small phase shift. It should be noted that the momentum 
equation including nonlinear term will be used for rest of the simulations in this 
paper. 

4.3. Effects of the assumed velocity profiles 

The coefficients of the momentum equation are based on the velocity profile that 
we assume when deriving pulsed flow equations. These coefficients are determined 
by the velocity profile over the cross-section of tube. Here we compare the classi- 
cal one dimensional model with steady Hagen-Poiseuille parabolic velocity profile 
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(7 = 2), the plug- like velocity profile (7 = 9), and the transient Womersley profile 
model proposed in this study. The effects of the assumed velocity profiles on the 
pressure waveforms are illustrated in Fig 12 at three different segments: Ascending 
Aorta, Abdominal Aorta I, and Right Femoral Artery. It can be easily seen that 
the mathematical model with plug-like velocity profile overestimates the pressure 
waveform due to the higher viscous drag along the whole arterial tree. It should 
be noted that all the other modeling options such as bifurcation point treatment, 
boundary conditions, and grid size are chosen identical to emphasize the true effect 
of the assumed velocity profiles. In order to see the effect of the assumed velocity 
profiles along the whole arterial tree, in Fig. |13| we perform corresponding surface 
pressure waveforms plot along the arteries form the Ascending Aorta to the Right 
Posterior Tibial Artery, colored blue in the Fig. [5] It is clear that choosing of the 
velocity profile makes a large difference in arterial system modeling, and we present 
a global model with smoothly changing velocity profile along the arterial tree. How- 
ever, if one would like to choose one of the Hagen-Poiseuille profiles, we conclude 
that the parabolic velocity profile seems a better choice to capture the waveform 
correctly instead of plug-like velocity profiles. Since the (3 parameter are less than 
22 in majority part of the arterial tree network, choosing plug-like velocity profiles 
(giving /? = 22) overestimates the pressure waveforms through the peripheral parts 
arterial tree. 
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Fig. 11. Effect of the nonlinear term in the momentum equation on the flow waveforms. Linear 
model assumes coefficient of the nonlinear term a = 0. 



180 r 




Fig. 12. Computed pressure waveforms with Hagen-Poiseuille and Womersley velocity profiles at 
the Ascending Aorta (Asc A), Abdominal Aorta I (Abd AI) and Right Femoral Artery (RFA). 
HP(7 = 2) assumes the parabolic velocity profiles corresponding a = 4/3, and /? = 8; HP(7 = 9) 
assumes the plug-like velocity profiles corresponding a = 11/10, and /3 = 22; and Wo assumes 
the functional coefficients (changing along the vessels in whole network) given by Eq.[9]for a and 
Eq.[l0]for /3. 
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Fig. 13. Surface plots of the pressure waveforms along the Ascending Aorta up to the Right 
Posterior Tibial Artery (including 14 segments); (a) Hagen-Poiseuille parabolic velocity profile 
(/3 = 8), (b) Hagen-Poiseuille plug-like velocity profile (/3 = 22), (c) Womersley profile (/3 from 
Eq. \10\ by varying /3 values through the vessels in the network. 

4.4. Effects of the bifurcation point modeling 

In previous studies, the boundary conditions at the bifurcating points in the network 
are mostly concentrated on the conservation of mass flow rate at that point along 
with the conservation of pressure or total pressure. Additional equations to close 
the nonlinear system of equations at those points come from either Riemann invari- 
ants or discrete continuity equations as explained before. Therefore we compared 
these four strategies to quantify the relative modeling difference. The computed 
pressure waveforms for Abdominal Aorta I and Right Femoral Arteries are shown 
in Fig. [14] Similarly, computed flow waveforms are also illustrated in Fig. [l5j The 
results demonstrate that the modeling with pressure or total pressure strategies at 
the bifurcation points effects the waveforms considerably. However, the complemen- 
tary equations coming from the Riemann invariants or discrete continuity equations 
does not effect the results. 

4.5. Effects of the outflow boundary conditions 

The imposed distal boundary conditions depend on the type of lumped parameter 
model used at the terminal points. Various lumped parameter outflow boundary 
condition models including Eq. [28] and Eq. [29] are available for arterial tree sim- 
ulations^ 30 * 38 * 46 *. When we analyzing the effects of certain modeling parameters in 
previous sections, we have used the pure resistance model. Here, we compared the 
three-element Windkcsscl and pure resistance models for the terminal points in or- 
der to see the compliance effects beyond the terminal points. The physiological data 
for terminal segments, including the Rt and Ct values, was tabulated id 34 ! 37 ! 41 !, 
and the constant value of R\/Rt = 0.2 was used for Windkessel outflow bound- 
ary conditions. Therefore, the biggest part of total resistance (80%) is attributed 
to the effect of the small vessel at the distal microvasculature. The compliances of 
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distal arteries (16.5%) were estimated to give a total arterial compliance of about 
1.0 mL/mmHg. In order to show the effects of these two models, the pressure and 
flow waveforms are compared at three locations. As shown in Fig. [16] and Fig. [17} 
the treatment of outflow boundary conditions have significant effects for modeling 
of arterial tree. 



5. Conclusion 

An improved model has been developed for the one-dimensional theory of flows 
through elastic tubes based on transient Womersley profiles changing continuously 
via the shape parameter. The model equations, the Pulsed Flow Equations (PFE), 
are a set of coupled partial differential equations that capture features particu- 
larly relevant to internal flows through flexible elastic tubes, such as flows in res- 
piratory and circulatory systems. The equations are an extension of the standard 
one-dimensional fluid flow equations that are able to capture some additional two- 
dimensional transport, viscous, and branching effects. The improvement of the one- 
dimensional theory by introducing functional coefficients to the momentum equation 
suggests that the formalism provided here can be used as a more appropriate global 
macroscopic reduced-order distributed mathematical model for the pulsatile fluid 
dynamics in physiological systems. The PFE are applied to a model of the human 
arterial tree containing the largest fifty-five arteries to investigate the blood flow 
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Fig. 15. Comparison of the bifurcation point treatments summarized in Table 2 by showing the 
effects on flow waveforms. 



characteristics for various velocity profiles and bifurcation modeling strategies, as 
well as the different lumped parameter outflow boundary conditions. The equations 
are discretized by the Lax-Wendroff scheme which is second-order accurate in time 
and space. The results demonstrate that the effect of the assumed velocity profile 
are significant for arterial tree modeling. The results also show that both the bi- 
furcation point modeling and the treatment of the terminal boundary conditions 
have also significant effects on the results. The quadratic nonlinear velocity term 
in the momentum equation also has a minor importance in modeling. Finally, the 
formalism is not restricted to the human cardiovascular system. It can also be used 
to compute other internal flows in flexible channels, including other physiological 
flows. Flows in the respiratory system, for example, and can be used as an efficient 
tool for computing the reduced-order macroscopic level of description of geometric 
multiscale models. 
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Fig. 16. Comparison of the lumped parameter terminal boundary conditions by showing the 
effects on pressure waveforms. 
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Appendix 



Since the contribution of the viscous dissipation term, 



, d 2 u 
dx 2 ' 



to the momentum 



equation is negligible, the governing PFE's given in Eq. [42] can be considered as 
a coupled hyperbolic system. The characteristics analysis of this system can also 
be performed by writing the system in non-conservative form. When Aq and cq are 
constant and a = 1, we can approximate this hyperbolic system in non-conservative 
form as: 



dQ +H dQ 
dt dx 



S 



where 



H = 



dF 
dQ 



A 



(46) 



(47) 



The eigenvalues of H are the roots of the characteristic equation \H — XI\. There are 
two distinct real eigenvalues A + = -0 + co(^ 1 ) 1 / 4 , and A~ = u — co(^) 1 / 4 , and hence 
the system is totally hyperbolic. The left eigenvectors are non-trivial solutions of 
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Fig. 17. Comparison of the lumped parameter terminal boundary conditions by showing the 
effects on flow waveforms. 



the system: 

[h h) 

Hence, the solutions are 



u - X A 

4^ A- 3 / 2 u- A 



lf=c Al /4 A- 5 /Ht 



[00] 



(48) 



(49) 



where is arbitrary. To obtain the Riemann invariants, q, we consider the following 
definitions: 



dq ± 
~dA 



' du ~ 2 

The integrating factor, u, is determined from the consistency condition: 

d(jj±lt) _ d(^lf) 



du 



dA 



(50) 



(51) 



By choosing l 2 = 1, and fx = 1, Riemann invariants which obtained by integrating 
the two equations given in Eq. ( 50 ) implying the initial conditions of u = at 
A = Aq, are given as: 



,±4c [l-(^oM) 1 / 4 ] 



(52) 
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Finally, the Riemann invariant form of the non-conservative hyperbolic system is 
given as 

^ (« ± 4c [1 - {A Q /A) 1/4 } ) + \ ± ^-(u± 4c [1 - {Ml A) x / 4 ] ) = - - ^ - (53) 



When the right-hand-side of the Eq. (53) goes zero, dq — on the characteristic 



line : 4f = A~, and dg + = on the characteristic line C + : ^j? = A + . Therefore, 
we can use these approximate Riemann invariants when we construct the boundary 
conditions. 
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